#include "kernel/chemsys.hpp"
#include "yocto/code/utils.hpp"
#include "yocto/math/kernel/algebra.hpp"


void ChemSys:: solveGamma( Real t )
{
	assert(ftol>0);
	bool converged = false;
	do
	{
		//----------------------------------------------------------------------
		//-- compute Xi (Chemical Jacobian)
		//----------------------------------------------------------------------
		computeGammaAndPhi(t,false);
		computeXi();
		
		//----------------------------------------------------------------------
		//-- compute Newton step
		//----------------------------------------------------------------------
		algebra<Real>::mul(dX,Xi,Gamma);
		
		//----------------------------------------------------------------------
		//-- update X
		//----------------------------------------------------------------------
		//std::cerr << "dX=" << dX << std::endl;
		
		//----------------------------------------------------------------------
		//-- check step
		//----------------------------------------------------------------------
#if 0
		bool full_step = true;
	CHECK_STEP:
		for( size_t i=1; i <= M; ++i )
		{
			assert( X[i] >= 0 );
			const Real dx = dX[i];
			if( dx > 0 && dx > X[i] ) // redundant...
			{
				for( size_t j=1; j <= M; ++j )
				{
					dX[j] *= 0.5;
				}
				full_step = false;
				goto CHECK_STEP;
			}
		}
#endif
		
		//----------------------------------------------------------------------
		//-- check convergence
		//----------------------------------------------------------------------
		converged = true; //full_step;
		for( size_t i=1; i <= M; ++i )
		{
			const Real dx   = dX[i];
			const Real xold = X[i];
			const Real xnew = xold-dx; //max_of<Real>(xold - dx,0);
			const bool ok   = Fabs( dx ) <= ftol * ( Fabs(xold) + Fabs(xnew) );
			X[i] = xnew;
			if( !ok ) converged = false;
		}
		
	}
	while( !converged );
	//std::cerr << "Xok=" << X << std::endl;
}
